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Abstract The paper considers random motion of a point on the surface of a sphere, in the ease where 
the angular velocity is determined by an Ornstein-Uhlenbeck process. The solution is fully characterized 
by only one dimensionless number, the persistence angle, which is the typical angle of rotation of the 
object during the correlation time of the angular velocity. 

We first show that the two-dimensional case is exactly solvable. When the persistence angle is large, 
a series for the correlation function has the surprising property that its sum varies much more slowly 
^ , than any of its individual terms. 

I ' In three dimensions we obtain asymptotic forms for the correlation function, in the limits where 

, the persistence angle is very small and very large. The latter case exhibits a complicated transient, 

O ■ followed by a much slower exponential decay. The decay rate is determined by the solution of a radial 

' Schrodinger equation in which the angular momentum quantum number takes an irrational value, 

namely j = ^{y/Tl - 1). 

' Possible applications of the model to objects tumbling in a turbulent environment are discussed. 



Keywords Diffusion, Ornstein-Uhlenbeck process 
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On , 1 Introduction 

If^ \ There are many contexts in which random motion is confined to the surface of a sphere. Examples 

■ include the motion of a unit vector n(t) indicating the orientation of an object tumbling in a turbulent 

fiuid flow, or the advection of a tracer in a thin, turbulent planetary atmosphere. In some applications 
it is sufficient to model the motion as diffusion on the surface of a sphere, which can be solved by 
noting that the eigenfunctions of the diffusion operator are spherical harmonics. In other applications, 
however, the angular velocity varies smoothly as a function of time, and the diffusive approximation 
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is not valid. Our objective is to obtain insight into random motion of a unit vector, by studying what 
is perhaps the simplest model. 

Smooth random motion on a sphere can be characterised by a dimensionless parameter which 
measures the typical angular distance through which the point has turned in the timescale for relaxation 
of fluctuations of the angular velocity. Due to an analogy with the concept of persistence length in 
polymer physics [T], we term this parameter the persistence angle: it will be denoted by (3. It is 
desirable to have a model which is a physically well motivated description of some processes, in which 
/? appears as a parameter. This paper discusses such a model, which is an extension of the well-known 
Ornstein-Uhlenbcck process [H[3] to describe motion on a circle or a sphere. 

We start by describing the standard Ornstein-Uhlenbeck equation on a line, in order to intro- 
duce some notation and elementary ideas. The Ornstein-Uhlenbcck process is a stochastic differential 
equation for the time-dependence of a variable v{t): 

V = -jv + V2Dr]{t) (1) 

where r]{t) is a white-noise signal, satisfying 

(rKt))=0, {iimt')}^S{t~t') . (2) 

Throughout this paper, (X) is the expectation value of X. The process equilibrates to a statistically 
stationary state, characterised by the following correlation function: 

{v{t + At)v{t)) = — cxp{^j\At\) . (3) 
7 

Equation ([T]) may be considered purely as a model for the fluctuations of velocity u of a particle, or 
it may be combined with the equation x = v to give a model for the displacement x. The motion 
in space is ballistic when viewed on short timescales, but on long timcscales it is diffusive, with the 
displacement Ax satisfying (Ax) = and {Ax'^) ~ 21)1, with spatial diffusion constant V: 

V^\r dt {v{t)vm = § . (4) 

^ J -oo 1 

To generalize the Ornstein-Uhlenbeck model on a line to a circle, we simply replace the variable v by 
the angular velocity, u, and the displacement along the line, x, by the angle along the circle, d. In two 
dimensions, we are considering the Ornstein-Uhlenbeck process on a manifold which has the simplest 
closed topology, namely a circle. 

If the equation of motion ([T|) is interpreted as a description of an angular velocity uj, then D has 
dimension [D] = T^''. The damping rate has dimension [7] = T~^. Thus, the problem is completely 
characterized by one single dimensionless parameter constructed from D and 7; we take this to be 

-/I 

The parameter ji has the following simple physical interpretation. The typical angular velocity is, 
according to (O, D/^, and its fluctuations occur on a timescale 7^^. The typical angle of rotation 
over the correlation timescale of the angular velocity is then AO ^ ^JD/^/^ = /3, so that /3 does 
correspond to a persistence angle. 

In this paper we show how to compute statistics characterising both circular and spherical Ornstain- 
Uhlcnbeck processes, such as the correlation function 

C{t) = (n(t) . n(0)) . (6) 

The solution to this problem has very different properties in two and three dimensions. 

In two dimensions, exact formulae are obtained for C{t) by expressing this correlation function in 
terms of the eigenvalues and eigenfunctions of the Fokker-Planck equation describing the probability 
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density of the angular velocity and angle variables. We find that the correlation function C{t) can be 
expressed as a series: 

oo 

C{t) = exp(-/3^t) ^ C^(/3) exp(-7V7t) (7) 

where the coefficients Cm{P) are determined explicitly, and satisfy Cjv(O) = 5^,0- We show that 
series given by Eql7] can be summed exactly, which reduces the correlation function to a very simple 
analytic form. There is, however, a surprising feature of this series expansion which deserves comment. 
Although the series ([7]) is convergent for all t and for all /3, in the limit as /3 — )- oo the structure of ([7]) 
appears to be hard to reconcile with physical expectations. We expect that when /3 ^ 1, the object 
rotates with almost constant angular velocity, and the decay of correlations results from motions with 
different angular velocity getting out of phase. This occurs on a timescale yJ^/D = 1/(7/?), iarger than 
the timescale for decay of the leading factor in ([7]), which is 1/(7/?^). The implication is that, when 
/3 3> 1, the summation over the exponentially decreasing terms approaches an increasing exponential, 
which almost cancels the rapid decay of cxp(— /3^7t). We show how this behaviour is realised, by 
taking very large coefficients with alternating signs. This effect is analogous to a phenomenon known 
as superoscillation where a Fourier sum oscillates faster than appears to be possible due to its 
bandwidth. In the case where /3 <i; 1, the motion of n{t) is well approximated by diffusion on the 
circle, and in that limit ([7]) yields C{t) = exp(— _Dt/7^), as expected. 

The extension of the problem to three-dimensions requires the introduction of a vector to describe 
the angular velocity. In addition, the closed manifold that describes the configuration of the system 
has a significantly more complicated topology. It turns out that in three-dimensions, the Ornstein- 
Uhlenbeck model of dynamic on a sphere does not appear to allow an exact solution, and we show 
that the spectrum has a very different structure from that of the two-dimensional case. When (3^0, 
the dynamics is easily understood in terms of diffusion on the surface of a sphere, for which the 
eigenfunctions are spherical harmonics. The solution has a much more complex behaviour in the limit 
as /3 — 00. As for the two-dimensional case, there is a transient which occurs on a timescale 1/7/?, but 
this transient is followed by a slow exponential decay of the correlation function on a timescale ^--^1/7. 
We show that in the limit as /3 — >■ 00 the slowest-decaying mode contributing to the correlation function 
is obtained from a radial Schrodinger equation for which the quantum number j takes an irrational 
value. Physically, the two time-scale solution we find for the correlation function can be explained by 
noting that when the angular velocity is constant, the component of n parallel to w is also a constant. 
The existence of this invariance explains the slow decorrelation of n at long times. 

Because this problem is far from straightforward we start by considering the two-dimensional case, 
where the Ornstein-Uhlcnbeck process has a circular coordinate. In section [5] the correlation function 
^ is determined in closed form for the circular Ornstein-Uhlenbeck process. The spherical Ornstein- 
Uhlcnbeck process is discussed in section [31 where we show that the Fokkcr-Planck operator can be 
transformed into a quantum Hamiltonian for a spin coupled to a spherical harmonic oscillator. Section|3] 
also considers the symmtery properties of the Fokker-Planck operator and the use of the eigenfunctions 
of the three-dimensional harmonic oscillator as a convenient basis set. The asymptotic behaviour of 
the solution for both large and small /? is considered in sectional Section [5] contains some discussion 
of possible areas of application of the model. Technical details of matrix elements required in section |4] 
are discussed in appendix A. Appendix B discusses an alternative model for continuous random motion 
on a circle, which is required to support the discussion in section [S] 



2 Two-dimensional case 

2.1 Formulation and general solution 

We consider motion on a circle, with angular coordinate 0, angular velocity w, replacing equation ([Ij 
by: 

9 = Lu , io = --fuj + \/2DT]{t) (8) 

where r]{t) is a standard white noise signal, with statistics satisfying and where 7, D are the 
damping and diffusion constants. Our aim is to be able to compute correlation functions such as 
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(n(t) • n(0)), where n(<) is a unit vector with direction 6. This correlation function is obtained by 
computing {cos9(t)) subject to the initial condition 9(0) = 0, with w having its equilibrium distribution. 

Statistics of the model ^ are obtained by computing the joint probability density of 9 and uj, 
P{9,ijj,t), with an appropriate initial condition. This probability density satisfies a Fokker- Planck 
equation [5] 

dP d d d'^P - 

This equation is solved by determining the spectrum and eigcnfunctions of the Fokker-Planck operator 
F, satisfying F 'Pnm{9,'^) = ^nm^nmi9 , uj) . This operator is separable, and its eigcnfunctions are of 
the form ^n„i{9,uj) = exp(in6')?/'„m(a;). The functions ^pnmiuj) satisfy 

d 

7^[wV'nm(l^)] + D-^lprimi!^) ~ inUJ1pnm{uj) = A„mVnm(l^) • (10) 

To solve this equation, define the operators 

Fa^jd^x + Ddl, F{a)^Fo-ax. (11) 

The operator F(a) is not Hermitian. It is convenient to introduce the following transformation to a 
Hermitian form: 

2 

H{a) = exp(7xV4-D)F(a) exp(-7a;V4L') = Ddl - Jj^x^ + ^-ax . (12) 

Note that Hq = -ff (0) is an inverted harmonic oscillator, with eigenvalues —7™, m = 0, 1, .... It will 
be convenient to use the Dirac notation for functions which are acted on by linear operators, replacing 
the usual angular 'bra-kets' with rounded ones to avoid confusion with our notation for expectation 
values. Thus the eigenfunctions of Hq are denoted by vectors \ fm)- 

Ho \(p,n) = -im |(p,„) . (13) 

The eigenfunctions will be assumed to be normalised according to the quantum mechanical conven- 
tion, so that the integral of their modulus squared over all space is equal to unity, for example the 
vector |(/3o) is a symbolic representation of the normalised cigenfunction (7/27rD)^/'* cxp(— 7x^/41)). 
These eigenfunctions may be generated from \ipo) by repeated application of creation and annihilation 
operators, a"*" and a respectively. These operators are 

a- = ^(^.-5.) , « = + . (14) 

Note that a+ is the Hermitian conjugate of d. The annihilation and creation operators, a, a+ satisfy 

Hn = -ia^a, [H(3,a^]^ --fo^ , [Ho,a]=ia, [a,d+] = l. (15) 

These relations imply that if Ho(p{x) = A(/?(x), then a'^ip{x) and aip{x) are also eigenfunctions, with 
eigenvalues A — 7 and A + 7 respectively (with the exception of the ground state, which is destroyed 
by d). The following relations describe the action of d, d"*" on normalised eigenfunctions of Hq: 



a\^Pn) = Vn\(Pn-i) , a+|</J,i) = \/n + l\tpn+i) ■ (16) 
The eigenfunctions of Hq are orthogonal, because this operator is Hermitian: 

POO 

ifnlfm) = / dx (pn{x)if,n{x) = Snm ■ (17) 



The eigenfunctions and eigenvalues for other values of a can be obtained by considering a transforma- 
tion to a new coordinate y = x + xq, where the shift is xq = 2Da/j^. In terms of this new coordinate, 
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we find that H{a) is equivalent to Hq + Dc? j^'^. The eigenvalues Am(Q:) and eigenfunctions |(^m(a)) 
of /f(a) are therefore 

A™(a) = -7™ + —Y- . I'^m(a)) = T{2Da/-f'^) |(^™) (18) 

where T{X) is a translation operator, defined by T{X)f{x) = /(a; — X). This operator may be repre- 
sented as an exponential 

f{X) = exp(-Xa,) . (19) 
By expanding the exponential and comparing with the Taylor seriers, we see that this expression is 
consistent with the defining property that T{X) f{x) = f{x — X) for any function /(x). Comparing 
pUI and pT|) . we see that we require a — in, so that we shall require matrix elements of this operator 
for complex values of X. 

The eigenfunctions Ifm) of Hq correspond to eigenfunctions ^pm{x) of F{a), where: 

^™(x) = expi^jx^AD) f{2Da/j^) (p™(x) . (20) 

If a function f{x) is written as a linear combination of these eigenfunctions 

00 

f{x) = ^ ajn1p,n{x) (21) 

then, in order to use the orthogonality property (|17L we must multiply f{x) by cxp(7a;^/4D) and then 
act on it with r(— 2Da/7^), to obtain 

dx ifmix) f{~2Dah'^) cMlx^/^D) f{x) . (22) 

In order to evaluate these coefficients, it will also be useful to have an expression for matrix elements of 
the translation operator (Franck-Condon factors) in the basis of the harmonic oscillator eigenfunctions: 



Inm{X) = {Lpn\T{X)\Lpm) = dx Lpn{x)ipm{x - X) . (23) 

J —oo 

By a simple adaptation of an argument presented in [5], for n> m these can be shown to be given by 
where L^"* (x) is the associated Laguerrc polynomial [B] : 

k=0 ■ ^ ^ 

The case n < m is obtained by using the fact that Imn{X) = Inm{—X). 

Using these results and noting that comparison of (|T0|) and (|TT1) implies a = in, we see that the 
eigenfunctions of the Fokker-Planck operator F are 

V'„™(0, cj) = exp(in0) cM-W/^D)f{2inD/-i^)^^{i^) . (26) 

These eigenfunctions are complex-valued for n 7^ 0, but note that the degenerate eigenfunctions ipn,m 
and tp-n.rn are complex conjugates, so that we can form real- valued solutions of 0. A general solution 
of Q can be written in the form 

00 00 

n— — 00 m— 

where the coefficients a^m are determined by the initial conditions: using (j22p gives 

anm - 7^ / de exp{-m9) / du ^Mf{~2inDh^) exp{-^uj^ /AD) P{e,Lo,Q) . (28) 

Jo J -00 



6 



2.2 Evaluation of correlation functions 

In order to evaluate the correlation function C{t) defined by ([§]), we consider the distribution of 6 with 
the initial condition that the initial orientation = 0, say) is known, but the angular momentum 
distribution initially in equilibrium: 

Pi0, ^' 0) = y ^ eM-l^Vm S{0) = (^) S{0) eM-l^Vm VoH ■ (29) 
Using p8[) . the coefficients in ((27| are 



1 / T \ 1/4 Z""" ^ 1 / T \ 1/4 

(^) J du:^MT{-2inD/-f'-)M^) = ^{:^) Irnoi-2inD/-f') . (30) 



2n \2nD 
The correlation function is 



C{t) = (n(t) • n(0)) = (cos 61) 

d9 cos 61 / do; P{9,uj,t) 

J — oo 

= 27rcxp(-L»t/72) ^ Re oi™ exp(-777ii) / dw cxp{-joj'^/4:D)f{2iD/'-f'^)if,n{oj) 

m=0 L -^-oc 

oo 

= exp{-Dt/j^) IrM{-2\Dh^) In,ni2iD/j^) exp{--fmt) . (31) 

Tn=0 

Note that ([21]) implies that the only associated Laguerre polynomials which are required are Lq"-* {x) = 
1, so that the Franck-Condon factors in (PT|) are 

/™o(2ii?/7') = \A(-i/3)'" exp(/3V2) = /o™(-2iD/72) (32) 
V m\ 



where /? = y^D/j^ is the dimensionless parameter defined in ([5]), so that the correlation function in 
(ED is 



C{t) = cxp{-Dt/j^) exp(/32) ^ ^AL exp(~m70 . (33) 

^-^ ml 

m- 

Using the identity 



2\m 



°° ( — x)™ 

; — exp(— am) — exp[— a; exp(— a)] (34) 



m 

m=0 



the correlation function can be obtained is closed form: 



C{t) - c(7t, /W) , c(t, /3) = exp[/32(l - r - exp(-T))] . (35) 

Figure [T] shows a comparison between the correlation function obtained by numerical averaging for 
the circular Ornstein-Uhlenbeck process, for three different values of (3. In each case the results are 
compared with the theoretical expression, equation (|35p. and the agreement is excellent. 
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2-dimensions; |3 = 1/12 



2-dimensions; p = 1/12 




5, 0.1 



0.01 




A 0.1 
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0.001 



0.1 



5, 0.01 



0.001 



0.0001 



2-dimensioiis; p = 3/4 




2-dimensions; p = 12 




Fig. 1 Correlation function for the circular Ornstein-Uhlenbeck process, for three values of — D/j'^: 
P"^ = 1/12 (upper row), — 3/4 (middle row) and — 12 (lower row). The simulations (full curves) show 
excellent agreement with the theoretical results, equation (f35)l. shown as a dashed line. 



2.3 Discussion 

In the limits /3 — > and /3 — >■ oo the correlation function approaches limiting forms which are expo- 
nential and Gaussian, respectively: from (|55|) we find 



C{t) - ey.Y'i-Dtl-f'^) , /3 < 1 
C{t) - ey.Y>{-De/2-i) , /3 > 1 



(36) 
(37) 



It is instructive to consider how these limiting cases arise. In the case where ^JdJ^ <C 1, the partial 
probability density P(0, £) satisfies a diffusion equation 



dP d'^P 
= V 

dt de^ 



(38) 



8 



which is to be solved with the initial condition P{9,Q) ~ 5 {9). The solution of this equation on the 
circle is 

oo 

P{9,t)= ^ a„ cxp(iTO6') cxp(-m22?t) . (39) 

m— — oo 

The initial condition gives = l/27r, so that 

^ oo .271- 

(cos6'(t)) = — V cxp(-7n^Pt) / d6' cos(6')cxp(im6') 

m=-oo 'J 

= exp(-2?t) = cx^{-Dt/-i^) (40) 



in agreement with (p6)) . In the opposite limit, ^jDj^ ^ 1, the particle rotates around the circle at a 
rate which is equal to its initial angular velocity. The equilibrium distribution of angular momentum 
has density 

^^""^ " cxp(-c^^/2i?) . (41) 

After time t the angle is = wt, so that the probability distribution of the angle is 



P{(^. t) = \j eM^O'l/^Dt') • (42) 

The correlation function is then 

{cos9{t)) = y dfl cos(0)exp(-02^/2i?^2) 

= exp(-£)tV27) (43) 

in agreement with (|37p . 

The introduction mentioned that there is a surprising aspect to the behaviour of the series (|33p in 
the limit as /3 — >■ oo. Despite the fact that every term in this series has an exponential decay with a 
timescale shorter than l/j/3'^, the exact evaluation of the sum of this series, equation (|55|). decays on 
a much longer timescale, 1/(7/?). 



3 Three dimensional case 

3.1 Formulation and Fokker-Planck equation 

We consider a unit vector n{t) evolving according to the equation: 

— =uj An (44) 

where u) is the angular velocity vector. The components uji of the angular velocity are determined by 
independent Ornstein-Uhlenbeck equations: 

= -^cj, + V2Drj,{t) (45) 

at 

where the i]i{t) are independent white noise signals, with statistics specified by equation This 
process is described by a Fokker-Planck equation for the joint probability density of n and u>. The 
general form for the Fokker-Planck equation in a space with coordinates Xi is [3] : 

2—1 t—1 J — 1 
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The velocities Vi and diffusion coefficients Dij are defined in terms of tlie expectation values of incre- 
ments Sxi in time St by writing (Sxi) ~ ViSt and (SxiSxj) ~ 2Dij5t. In our case (xi, X2, . . . , xg) = x ^ 
(n,u)) = (ni, n2, 71-3, wi, a;2, wa). The velocities are 

Vi = eijki^jUk , i = 1, 2, 3 , Ui = -7CJi , i = 4, 5, 6 (47) 

and the diffusion coefficients are 

D,, = DS^j , i = 4, 5, 6 , = , i = 1, 2, 3 . (48) 

The Fokker-Planck equation is therefore 

f - - ^.^-^(-.-^) + 7^(^.P) + (49) 

where in equation (j49p . as well as in the equations below, repeated indices are summed over the values 
1, 2, 3. The Fokker-Planck operator can also be expressed in the form 

T ^ "fdiUJi + Ddidi ~ oJiJi (50) 

where Ji are components of an angular momentum operator, defined by 

J^ = 'vkn,^^ (51) 

and where di = d/duii. Note that this definition differs from that which is commonly used in quantum 
mechaincs texts (such as [7]) by a factor of i = \/— T. The Fokker-Planck operator ([50)) can also be 
expressed in the form: 

J- = - u; • J , J'o = 7V • a; + L>V • V . (52) 

where V is the gradient in the angular momentum space. 
The variables may be made dimensionless by using: 

< = 7^ (53) 

—uj, . (54) 
2D ^ ' 

In dimensionless form the Fokker-Planck equation reads: 

dtP = d,{u,P) + ^/2^3e,,k^k^n, {n,P) + ]^djP (55) 

where the variables in equation (|55p are dimensionless, as defined in equations (|53l54p : for simplicity, 
the overbars have been omitted. 



3.2 Symmetry analysis 

Consider the symmetry properties of the Ornstein-Uhlenbeck operator describing the evolution of 
the joint probability density function P{n,<jj,t). The physical properties of the system are invariant 
under the rotation group, in the sense that an arbitrary rotation of both n and ui leaves the problem 
unchanged. We define angular momentum operators which generate these rotations: 

Li = €ijkiOjdi^^ 

Ji = eijkTijdn^ (56) 

(the operators Ji were already considered in ([FT]) ). These operators are related to the operators used 
in quantum mechanics by a factor i, that is L — iLqm where Lqm is the angular momentum operator 
defined in standard texts such as [7]. Thus, adapting standard results [7], the commutation relations 
of the operators are: 

[Li,Lj] — —€ijkLk (57) 
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and similar relations for Ji. It is also useful to note that: 

[Li.ujj] = -eijkUJk ■ (58) 

The term a; • J is invariant when rotating simultaneously n and a;, and using (j57p and ()58[) it is easy 
to check that: 

[J + L,uj-J]=0. (59) 



[y,uj-J]=0. (60) 



In the same way, it is easy to see that 



It is clear that these are also symmetries of the other elements of the Fokker-Planck operator in ([55|) . 
Because the system is invariant under rotation simultaneously of a; and n, the Fokker-Planck evolution 
operator commutes with the total angular momentum, as well as with the magnitude of the angular 
momentum of n{t). The set of conserved quantities is therefore 

k = L + J, f (61) 

where L acts on the u) variables and J acts on the n variables. The corresponding quantum numbers 
are denoted j, mj for J, I, nii for L, and fc, rrik for K. 

Thus, the Fokker-Planck operator has a block-diagonal structure, where the blocks are labelled by 

a given set of values of j, k and nik, where -I- 1) is the eigenvalue of J , mk the eigenvalue of K^^ 

- 2 

and ~k(k + 1) the eigenvalue of K . These symmetries indicate that spherical polar coordinates and 
spherical harmonic functions will prove useful. We use {9,(j)) as the polar angles for n and {uj,9',(j>') 
as the spherical polar coordinates of u:. The spherical harmonic functions Yim{0, (j)) will be denoted by 
Dirac state vectors, writing |Yi„i) when the arguments are {0,(j)) and |i^/„) when the arguments are 

Further constraints due to symmetry considerations can follow from the initial conditions. For ex- 
ample, the initial condition for evaluation of correlation functions has the variable lj in the equilibrium 
state, which is spherically symmetric with / = (and which is easily seen to be a Gaussian function of 
w). The value of n is set equal to one particular vector, say e^. The initial condition is therefore 

P(n,a;,0) = V4^(^)^^'exp(-7c^V2^) >oo(e',</'') ^(n - e,) . (62) 

The (5-function distribution can be resolved into a sum the spherical harmonics, with j = 0,1,2,... 
but with nij = in each case: 

<5(n-e,)=^-^P,(n.e,) (63) 

where the functions Pj{x) are Legendre polynomials. This means that we can confine attention to the 
TUk = subspace when evaluating an equilibrium correlation function of n(i). 

Furthermore, in the case of the simplest correlation function, C(t) = {cos 9), the quanitity being 
averaged spans only the j = 1 subspace, so we can confine our attention to the j = 1 subspace. Thus in 
order to evaluate C{t) wc must consider the j = 1, nik = subspace. Furthermore, the initial condition 
has / = 0, so the total angular momentum is fc = 1. The values of and are constants of the 
motion, indicating that wc must consider only solutions with fc = 1 and j ~ 1. Setting k = 1 and 
j = 1, the triangle relation for L = K — J is \k~j\<l<k + j, that is, / = 0, 1, 2. 

Wc can construct functions of the angular variables with definite values of J^, i^, and Kz, 
labelled by quantum numbers j, I, k, mk- Let these functions be denoted by Tj^i^k,mk{9,4'T9' and 
we assume that these functions are normalised so that they form an orthonormal set, with the usual 
integration measure for a cartesian product of two spherical surfaces. The symmetry considerations 
discussed above imply that the solution in the j = 1 subspace may be written in the form 



Pi(n,u;,i) = 5^V'i(w,i)n,u.o(e,'^,^','/'') . (64) 

1=0 
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This result shows that, if we are concerned with evaluating the correlation function C(t) = (cosO), 
then symmetry considerations reduce the Fokker-Planck equation to a system of three coupled ordinary 
differential equations. More generally, the calculation of correlation functions such as {Pj (cos 9)) implies 
values of I in the range < I < 2j, so the solution requires 2j + 1 functions of cj. 



3.3 Harmonic oscillator basis 

It will also be useful to consider the exact solution of the Ornstcin-Uhlcnbeck process describing the 
evolution of the angular momentum independent of evolution of n. This will be related to the 
three-dimensional quantum spherical harmonic oscillator. 

The operator J-q in (j52p has a structure which is closely related to the harmonic oscillator of quan- 
tum mechanics. It is convenient to transform the Fokker-Planck operator into a three-dimensional 
isotropic harmonic oscillator. We consider the operator 

7i = exp(x/2).Fexp(-x/2), ^ ^ lM±A±^ , (65) 

Using the dimensionless variables defined by ([5^ and ([5^ . we can express Ti in the form 



d+ • d + V2PJ -uj] =no- V2P7J ■ (66) 



where the components of d''' and d are, respectively, creation and annilhilation operators for the i 
degree of freedom, using the dimensionless variables defined by equations ([5^ [M]) : 



a. 



+ 



1 , „ , 1 



^ {uji - di) , ai = -j= [uii + di) . (67) 



The eigenvalues of Tio are — 7(fci + k2 + k^), where hi = 0, 1, 2, Note that all of the eigenvalues 

except the ground state are degenerate. 

As well as being separable in Cartesian coordinates, Hq is also separable in spherical polar coordi- 
nates: 

no = 7 [^V^ - + ^] (68) 
where the Laplacian operator may be expressed as 

Vl^dl + -d^ + ^L' (69) 

and L is, up to a factor i, the usual angular momentum operator acting on Sj^; its eigenvalues are 

+ 1), I = 0, 1,2, The degenerate multiplets can, therefore, also be resolved as states which 

- 2 

are eigenfunctions of L and L^, labelled by quantum numbers n,l,m (where + 1) and m are 
eigenvalues of L and respectively). The ground state eigenfunction of lio = — 7d^ • d is 

Vooo{^) = ( cxp ( - Y ) (70) 

and the other eigenfunctions of T-Lo are of the form 

^nlm{0j) = V^YirniO',(j)')uJ^Pnl (w^) V?OOo('^) (71) 

where Pni{x) is a polynomial of degree n, proportional to the generalized Laguerre polynomial 
^(/+i/2) rpj^g eigenfunctions are normalised in the usual way: 



dx :E2+2'exp(-x2)P„,(2:')P„'i(x2) 
-/ dj/y'+i/2cxp(-y)P„,(j/)F„,,(y) . (72) 
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Using the orthonormality relation for the generahzed Laguerre polynomials: 

Jo n\ 
we deduce the relation between P„i and the generalized Laguerre polynomials: 

Pnlix) = Mnl Lil+i\x) , = 



2n! 



(n + ; + i)r(n + ; + |) 



(73) 



(74) 



3.4 Equations of motion for J = 1 modes 

Here we consider how to write an equation of motion for the projection onto the j = 1 modes, which 
contribute to the correlation function C(t). The harmonic oscillator basis which was introduced in 
section [3.31 will prove useful here. 

In section 13.21 we showed how symmetry considerations constrain the angular dependences of the 
solutions. The angle-dependent parts of the solution are constructed from functions with known values 
of J^, L"^, and K^^ with quantum numbers j,l,k, nik- These functions may be expressed in terms 
of tensor products of spherical harmonics, writing 

|3^j,i,fc,mJ = XI ('i''^i!^2,m2|j,/,fc,mfc)|Y(i,r„J ® |y4,„J (75) 

where |Yii,mi) represents the spherical harmonic Yi^^mi {6, 4') which is a function of polar angles repre- 
senting the direction of n, and jYJ'^ represents Yi^^,ri2{(^' ■4''): which is a function of the polar angles 
for u). The coefficients {Ii,mi;l2,rn2\j,l, k,mk) are termed Clebsch-Gordon coefficients (see, for ex- 
ample, [Hl[Z])j and this representation is useful because the spherical harmonics have well-known and 
convenient properties. 

In particular, determining the correlation function C{t) requires a solution involving just three 
functions ipi{Lo,t)^ multiplying the angular functions Ti i i^, with / = 0,1,2 (see equation (j64p ). From 
tabulations of Clebsch-Gordon coefficients we find: 

|ri,o.i,o) = l^io) ® iFo'o) 

= /fi^i.+i) ® \y(.-i) - vi'^i-i) ® i^i'+i) 

|ri,2,i,0) = /^Ii"2,l) ® in-l) - l^2,0) ® \Ylo) + ® |yi', + l) . (76) 

In order to express the equation of motion in the j = 1 subspacc, it is necessary to rewrite the operator 
a; • J so that its action upon the spherical harmonics is explicit. To this end, consider the angular 
momentum ladder operators: 

J+ = {X + iJy) , J- = {Jx - iJy) (77) 

It is straightforward to see that: 

[J„J±]=±iJ±- (78) 

Thus, applying the operator to the eigenstate of Jz with a quantum number m, namely the spher- 
ical harmonic 1^;,™), one finds Jz^±|^,m) = i(?7i ± l)>/±|^,m)- The operator J+ thus increases the 
azimuthal quantum number by one, whereas J_ decreases the azimuthal quantum number by one. For 
completeness, the prefactor, up to a phase, can be obtained by expressing J as: 

/ = J+J- + Jl - Uz ^J-J+ + Jl + iJ, (79) 
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which immediately leads to: 



J+\Yi^m-i) = Wil + m){l-m+l)\Yi^m) 
J-\Yim) = Wil + m)il-m+l)\Yi^m-i) 



(80) 



With these results in place, we can express a; • J in terms of the operators J±. Elementary algebra 
leads to: 

UJ J = -{Ux - iUJy)J+ + -{UJ^ + iUJy)j^ + . (81) 

If one notices further that: 



/ 87r 

{ujx - iujy) = +LuJ — Yi^_i(0', (f)') 
I 4tt 

wd — Yi,o(6'',(/)') 



then, equation ([5T|) leads to 

/Itt / 1 



UJ ■ J — U!\ 



3 vy2 



(82) 



(83) 



The formulation of equation (1831) is useful to understand how the operator a; • J couples modes to each 
other. In particular, with the help of equations (j83p and (|76p we can determine the matrix elements 

(84) 



(^i,i,i,o|'^ • J|Ti,j,i,o) = wA. 



where i,j G {0, 1, 2}. Only four of these matrix elements are non-zero. After a lengthy but mechanical 
calculation we find the following values for the non-zero matrix elements: 



^01 
Al2 



Alo 
A21 



(85) 



It is convenient to use the connection with the spherical harmonic oscillator, and to replace the functions 
in (El by 

0(w,i) = exp(-wV4)V'K^,i) (86) 

(here we use the dimcnsionlcss variables (|53p . ([M)) ). Substituting (IMl) . ([55)) into the Fokker-Planck 
equation, multiplying by |Ti.j^i.o), and integrating over the product of two spheres, we obtain three 
partial differential equations for the three components coupling to the j = I mode. These equations 
are: 



d_ 
di 
d_ 
di 
d_ 
di 



Ciioj,t) = LiCi(w,t) + V2^Lj [^ioCo(w,i) +-4i2C2(w,i)] 

C2i0J,t) ^ 12(2(10, t) + V2^A2lCliLO,t) 



(87) 



where 



^ 2 



0^ + ^a.-^(^-c.^ + 3 
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4 Asymptotic properties of the correlation function 

In three dimensions we are only able to determine the spectrum of % (defined by equation (j66p ) 
by analytical methods in the limits /? — )■ and /? — ?• oo. This section considers various asymptotic 
approximations for the correlation function. 



4.1 Short-time limit 

The correlation function C(t) = (n(i) • n(0)) can be calculated by using the initial condition n(0) = es 
and then computing {cos 9). For short times, the polar angle is approximated hy — (a; J + w| 
0{t^). In the short-time limit, therefore 

C{t) ^{l-\e^ + ...)=l-\{ul+ u^Dt' + . . . (89) 

Using equation ([3]), we have {uSj) = D/^. The leading order behaviour of the correlation function is 
therefore 

C{t) = 1 + Oit") ^1-P^P + 0{t^) (90) 

7 

where t is the dimensionless time defined by ((53)) . This result is valid for all /3. 



4.2 Difi'usive (/3 0) limit 

In the limit Df-f^ ^ 1, the angle of n diffuses with diffusion coefficient V ~ The solution of the 

diffusion equation on the surface of a sphere is expressed in terms of spherical harmonics YimiO, 4>), for 
which the eigenvalues of the Laplacian are + 

00 / 

P{e, 0, t) = ^ cxp[-l{l + l)Dt/j^]Yir„{0, ct>) . (91) 

The required coeflicients are obtained using orthogonality of spherical harmonics: aim = Ym{^z)^m[), 
so that aim = (5,no'5ii-\/3/47r. The expectation value of cos{9) therefore has a simple exponential decay: 

/ 47r 

— (Yio(t)) = exp(-2I?t/72) = cxp(-2/32t) . (92) 



This approximation does not have the correct limiting behaviour as t — >■ 0, which is given by (j90p . The 
following approximation to the numerically determined correlation function approaches (|92[) at /3 — >■ 
for all t and has the correct quadratic behaviour at i = 0: 



4.3 Large j3 limit 

In the limit D /^^ ^ 1, the short-time behaviour of the correlation function is determined by tumbling 
motion with fixed angular momemtum. Consider the solution of the equation of motion h = w A n in 
the case where oj — Loe^j is constant and where the initial direction is no — e^^. The solution is 

n(t) = a + bcoswi + csincji (94) 

where 

a=(e^-no)e^, b = ng - a , c = ° \b\ (95) 

\a A o\ 
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are three mutually orthogonal vectors. It foUows that 

n(<) • n(0) = |no • e^p + cosa;i[l - |no • 
Now integrate over the distribution of angular momentum to obtain: 



(96) 



(n(i).n(O)) 



( 7 
\2TrD) 



cos ujt { 1 — 



\2ttD) 



dw. 



dw,, 



3/2 



= 47r 

_ 1 
^ 3 



2iiD, 
( 1 



\2tiD 
2 

^3 



1 - 



doJz exp(-7a;^/2L)) 

/ doj d9 2TTLj'^sm0cxp{-ju^/2D)[cos 
Jq Jo 

du} oj^ exp(— 7ti;^/2L)) 
exp{-Dt^/2j) . 



1 - 



cosut 



sin^ 6 cos ujt\ 



3/2 



/o 



+ - cos LUt 



7 



(97) 



This shows that in the three-dimensional case, when /? = Djy'- ^ 1 the correlation function decays 
to (n(t) • n(0)) = 1/3 on a rapid timescale, t\ = (/37)^^, due to motions with different frequencies 
getting out of phase. 

On longer timescales the value of a; fluctuates, and the correlation function then decays to zero on 
a slower timescale describing the decay of correlations of a;. A precise understanding this limit requires 
us to carry out a more sophisticated analysis, as will be done in section [4.41 below. Before addressing 
this issue we consider how n(t) behaves when uj varies slowly. We show that the angle between n and 
a; is an adiabatic invariant of the dynamics of (|33]). This shows that the decay of correlations of n(i) is 
governed by the diffusion of the direction of a;, implying that the timescale for the decay of correlations 
of n(<) is 0(7~^) in the limit as /? -> oo. 

To show that the angle d between n(t) and u) is an adiabatic invariant, consider the time evolution 

of 

f = n ■ u) = Ld cos = bjz (98) 
where uj = and where the second equality defines z ~ cosO. From (j44[) . the time derivative of / is 



dt 



du) 

dF 



5n{t) 



du) 
'dt 



(99) 



where Snit) oscillatates on a timescale cu ^ about a mean value which is equal to zero when u} is 
constant. When evaluating the drift of /, we neglect this rapidly oscillating term, and write 



dt 



z 

—ijj 



dijj 
dT 



dL^2 



2uj dt 



Alternatively, from the definition f = zuj, we find 



dt 



dt 



dio 
'd^ 



diU^ 



dt 2uj dt 



(100) 



(101) 



Comparing (|100p and (jlOip we see that 2 = 0, implying that 9 is an adiabatic invariant, as stated 
above. 

We have argued that in the limit as /3 — >■ oo the direction of n{t) is determined by the evolution of 
a;. This indicates that we can determine the rate of decay of C{t) by determining the rate of decay of 
correlations of the angular momentum. The direction vector of the angular momentum, exhibits 
diffusion on the surface of the unit sphere with a diffusion coefficient which we determine shortly. 
In a diffusiuon process, the expectation value of the spherical ha rmoni c cos 9 decays exponentially: 
(cos0) ~ exp{—2Vi^t). We might therefore expect that when /3 = y/D/j^ 1, the correlation function 
is well approximated by 



{n{t) ■ n(0)) ^ - exp(-At) + - (^1 - ) exp{-Dty2j) 



(102) 
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and that A = 21)^. This argument is, however, not satisfactory. When the deviation from equilibrium 
of the distribution of n(t) is well approximated by the spherical harmonic cos 6*, the corresponding 
distribution of e^^ might be a quite different combination of spherical harmonics. This question is most 
effectively addressed by the application of group theory to the Fokker-Planck equation, using results 
from sections 13.21 and 13.41 above. 

We conclude this section by determining the diffusion coefficient for diffusion of the direction of u). 
This direction is the unit vector = lj/lu, which diffuses on the unit sphere with diffusion coefficient 
V^. The change in the e^^ in a short time St is 

6e^ = [ujSu — (a; • Suj)uj] (103) 

so that the expectation value of the square of the rotation angle is 

{Sel)^^^{5u.f)^^. (104) 

The required diffusion coefficient is obtained by averaging this over the known distribution of u;: 



ASt J uj^ 

n(^Y" rd.W -"-^f/^°' .., ,105) 



JO ^ 

The correlation function for the direction of the angular momentum therefore deays as 

(eUt) ■ eUO)) = exp(-22?„t) = exp(-27t) . (106) 

Numerical evidence indicates that the correlation function C{t) decays at a different rate when /3 ^ 1: 
the decay rate in (jl02p is found to be A ~ I.567, rather than A = 27. The difference is explained in 
the following section, 14.41 



4.4 Asymptotic solution of slowest mode 

Here we use results derived from symmetry considerations (in sections 13 . 21 and 13 .4|) to identify the slow- 
est decaying modes in the limit as /3 — > 00. The objective is to determine solutions of ([57]) which decay 
exponentially in time, so that Q{uj,t) = exp(— At) aj{ui). The functions 0^(0;) satisfy the eigenvalue 
equation 



The structure of the operator 



Li ai(w) + \/2i3ll> [Aioao{io) + yli2a2(w)] = A ai(aj) 
L2a2{uj) + V2PLuA2iai{uj) — A 02 (w) 
implies that 



hm^^ 



(107) 



(108) 



for some constants Cj . 

We are interested in the most slowly decaying solutions, which requires determining the eigenvalue 
A with the largest real part. The structure of (|107p suggests that the decay rates are expected to 
increase in proportion to /3 as /? — >■ 00. However, our discussion in section 14.31 indicates that there 
should be eigenfunctions which have a slow deacy rate, A = 0(7) as /3 — >■ 00. In order to identify these 
slow modes, note that the matrix Aij has a null eigenvector: this matrix is 



(109) 
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which has a null vector, (1,0,2). The only way to obtain an eigenfunction with an eigenvalue which 
remains bounded as /3 — > oo is to assume that throughout most of the range of w, we have 



^10 ao(w) + Ai2 02(1^) ~ . 
Using this approximation, the equation for 02 (w) can be written 



1/2 + (Oic 



02 = (1 + p) Xa2 



where 



However, we have 



so that a2{uj) is an eigenfunction of 



A21A12 ^ 1 
AoiAio ~ 2 ■ 



6 1 



(110) 

(111) 
(112) 
(113) 
(114) 



l + pu;2 

with eigenvalue A. This is an operator of the form (jll3p with an angular momentum quantum number 
j which satisfies j{j + 1) = 6/(1 + p) = 4. that is. 



Vvf- 1 



(115) 



Thus it is argued that, in the limit as /? — > 00, there exist modes for which the eigenvalues are 0(7), 
which satisfy a radial equation with an irrational value of the angular momentum, given by (|115p . It 
remains to identify the eigenvalues associated with this equation. The first step is to factor out the 
boundary condition at w = 0, writing 

a.j{uj)=uj^<P{uj) . (116) 

Then it is useful to remove the Gaussian factor from the solution, writing •P{uj) = exp(— w2/2)(/)(a;2). 
The equation for (/)(w) (with u = w^) is: 



u^" + {] + 3/2 - u)(j)' + -^^—^0 = 



(117) 



Well behaved polynomial solutions exist only when (A — j) = 2n {n integer). Inserting the factor of 
7 which is required when we return to dimensioned equations, we conclude that the cigenvlaue which 
gives the slowest rate of decay is 

7(\/T7- 1) 



A = -7J = 



and the full set of eigenvalues of the problem defined by equation (|117p is 

A„ = -7(.?+2n) , 71 = 0,1,2,... . 



(118) 
(119) 



The eigcnfunctions can also be expressed in terms of generalised Laguerre polynomials, so that the 
Fokker-Planck equation may be regarded as exactly solvable in the limit as /3 — >■ 00. 

We have imposed two apparently incompatible conditions on the solutions aj{uj), namely (jlOSp 
and (|110p . We conclude this section by considering how these are reconciled. Consider the nature of 
the solutions aj{Lo) close to w = 0. Dimensionally, the problem is analogous to one with the following 
structure: d^^/dw^ + /3a;<^ = 0. The derivative terms becomes dominant when a; < 13'-^^^, so that the 
approximation (|110p fails and ()108p becomes applicable when w/?^/^ ^ 1. The problem can be treated 
by applying standard asymptotic expansion methods |^ . Alternatively, one can determine the solutions 
of the problem by expanding the solution on a conveniently complete basis, and by truncation, reduce 
the solution of equations (jl07p to a matrix equation, as explained in Appendix A. 

Figure [5] shows the numerically computed correlation functions for the three-dimensional Ornstein- 
Uhlenbeck process, compared with various asymptotic approximations. We also investigated the spec- 
trum of the Fokker-Planck operator numerically. Using the eigcnfunctions of the spherical harmonic 
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Fig. 2 Correlation function for the spherical Ornstein-Uhlenbeck process, for three values of = -D/7'^: 
= 1/48 (upper row), = 4/3 (middle row) and — 192 (lower row). The data for/?^ = 1/48 show good 
agreement with the exponential approximation to the correlation function, equation H92[l . which is applicable 
in the li mit as /? — >■ 0. The linear plot for 0^ = 192 shows good agreement with the transient described by 
equation (|102|l . and the logarithmic plot demonstrates that the correlation function C(t) behaves at long times 
as C{t) oc exp(— (vTZ— l)^jt/2), as predicted by equation {TTS]). 

oscillator as a basis set, this operator can be represented by an infinite-dimensional matrix. The formu- 
lae for the matrix elements are given in appendix A. We find that the spectrum of finite dimensional 
truncations converge as the size of the basis set increases, and we identify the converged eigenvalues 
with elements of the spectrum of the Fokker-Planck operator. Figures [3] and |4] illustrate the dependence 
of the eigenvalues upon /3. 



5 Applications to random tumbling 

This paper has described a model for the statistics of a unit vector n(i) moving randomly but smoothly 
over the surface of a sphere. The model is a generalisation of the Ornstein-Uhlenbeck process to a 
spherical geometry. It has the merit of being susceptible to analytical treatments, including an exact 
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Fig. 3 The six slowest decaying eigenvales from the spectrum of the Fokker-Planck operator of the spherical 
Ornstein-Uhlenbeck process as a function of /3. Left; real part of the eigenvlaues, showing the separation 
into levels which diverge as /3 — >■ oo and levels which approach a finite limit. Centre; real parts at greater 
magnification. At very small values of /3, all eigenvalues are real. Pairs of branches collide at a finite value of /3, 
giving rise to pairs of complex congugates eigenvectors, whose real parts are shown by the dashed lines. Right: 
positive imaginary parts of the eigenvalues. 




Fig. 4 The r eal-valued eigenvalues approach the spectrum A„ = j + 2n in the limit as /3 — > oo, in accord with 
equation (|118|l . The right-hand panel shows the convergence towards this limit as /3 — s> oo. 



solution in two dimensions. It is of interest to consider the possible applications of this model, and the 
extent to which the model provides an accurate description of various systems. 

The spherical Ornstein-Uhlenbeck model could be used to describe the tumbling of an object in 
a turbulent fluid flow. Examples include the rotational motion of rocks or dust grains in turbulent 
circumstellar discs [l^, or of ice crystals in a convecting atmosphere [11]. The rotational motion of 
such bodies can influence their growth by aggregation or their evaporation by exposure to a source of 
radiant heat. Recently it has become possible to make detailed experimental studies of the orientation 
of a neutrally buoyant sphere in a turbulent fluid by matching images of an irregularly painted ball 
to photographs taken with the ball in a well defined orientation [T2l[T3] . The direction n(t) of one 
axis through the sphere can be followed as a function of time. The model could also be used for the 
fluctuations of direction vector n{t) of a rod-like body in a turbulent fluid. 

These remarks raise the question as to whether the spherical Ornstein-Uhlenbeck model will give 
an accurate description of the tumbling motion of a body. In the case of a small body in a turbulent 
flow, we argue below that the statistical properties of the velocity gradients of turbulence appear to 
make this a very good model. The orientation of a small body in a turbulent flow with velocity field 
u{r{t),t) responds to the gradients of the the velocity field, evaluated along the trajectory r{t) of the 
body (which can be assumed to be advcctcd with the fluid). The velocity gradients form a matrix A(t), 
with elements Aij{t) = dui/drj{r{t),t). It is convenient to write A = f2 + S, where f2, the vorticity 
tensor, is antisymmetric and where S, the strain-rate, is symmetric. The equation of motion for the 
direction vector n(t) of a microscopic ellipsoidal object in a fluid flow was obtained by Jcffery [T3]. It 
can be written in the form 

^ = f2{t)n + 4^ [S(t)n - (n • S(t)n)n] (120) 
di -I- 1 
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where a is the axis ratio of the elhpsoid. The same equation of motion apphes to general axisymmetric 
bodies, provided they are small compared to any characteristic lengthscale of the flow [TS], but the 
relation between a and the shape of the object is not known in general. In the case of a spherical 
particle or other object with a ~ 1, the equation of motion (|120p is of the same form as equation (|33]), 
if we interpret the vorticity ui{t) as being the antisymmetric tensor corresponding to the axial vector 
with elements related by uji = Cijkf^jk- 

Furthermore, the Lagrangian correlation function of the vorticity in turbulent flows has been investi- 
gated using simulations of turbulent flows by several authors P^ITTlfTS] . The elements have mean value 
equal to zero and appear to be statistically independent. It was found that the correlation function of 
each elements can be fitted quite accurately by an exponential function: 



The decay rate 7v is of the order of the inverse of the Kol mogo rov timescale tk of the turbulent flow, 
which is the shortest timescale of the fluid motion: tk = ^JvJE^ where v is the kinematic viscosity and 
where £ is the rate of dissipation per unit mass. Numerical evidence indicates that 7v ~ 1/(8. Stk) at 
large Reynolds numbers [T8|. The diffusion coefficient can be related to 7v by various kinematic 
constraints (discussed in [TB]), giving = 7v/I2t^. These considerations suggest that the spherical 
Ornstein-Uhlenbeck model should describe tumbling of a small object in a turbulent flow, with a 
'universal' value for the persistence angle 



However, we cannot conclude that correlation function C(t) = (n(i) • n(0)) for objects in a turbulent 
fluid will correspond to that of our spherical Ornstein-Uhlenbeck model, because vorticity may have 
very different temporal variation, compared to the Ornstein-Uhlenbeck model, but still have the same 
correlation function. This point is illustrated by a calculation in appendix B, where we analyse a model 
for random motion on a circle, in which the angular velocity is determined by a telegraph process 
(examples of this type of model are considered in [1^120) ). The telegraph model has a correlation 
function of angular velocity which is an exponential function, equivalent to that of the circular Ornstein- 
Uhlenbeck process. However we show that the correlation function (n(t) • n(0)) is very different for the 
two models. We conclude that the extent to which the spherical Ornstein-Uhlenbeck process is a good 
description of tumbling in a turbulent fluid must be tested by numerical simulations of turbulence. 
Our own numerical investigations on the tumbling of microscoipic particles in turbulence indicate that 
the spherical Ornstein-Uhlenbeck model is not a good model for their correlation function [TB]. This 
conclusion is consistent with studies by Shin and Koch [21], who presented data for (n(i) • n(0)) in 
simulations of rod-like objects in fully developed driven turbulent flows. 

It is only in cases where the statistics of the angular momentum fluctuations arc a precise match to 
the spherical Ornstein-Uhlenbeck process that reliable predictions can be made about the correlation 
function defined by ([H]). In the case of an object tumbling in a very dilute gas, such as a small rock 
in the circumstellar disc of the star, the spherical Ornstein-Uhlenbeck model may be a very good 
description of the evolution of the angular momemtum. Bombardment by microscopic dust grains can 
provide random impulses which change the angular momentum in the same way as the white-noise 
fluctuations in (|45p . And the damping due to motion in an extremely dilute gas is proportional to the 
relative velocity [121 1 consistent with the linear damping term in equation (|45p . We conclude that a 
small body tumbling in a very dilute gas is one example where the spherical Ornstcin Uhlenbeck model 
is an excellent description of a physical process. 



6 Concluding remarks 

Our study was motivated by recent works, which have characterized the orientation of particles trans- 
ported by a turbulent flow [T^[T51[Tl 11231124) . These processes are naturally described by the random 
motion on a spherical surface. The work here has focused on arguably the simplest model for smooth 



(%(i)%(i')) - — cxp(-7v|t-i'|) . 



(121) 




(122) 
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random motion on a sphere: the direction n rotates with an angular velocity u;, which evolves according 
to an Ornstcin-Uhlenbeck process. 

The solution of this model depends on only one dimensionless parameter, which we refer to here as 

the persistence angle, which characterizes the rotation occuring during one correlation time of uj. 
We have characterised evolution of n(t) by analysing its correlation function (n(O)n(t)). 

In two dimensions, we have obtained an explicit expression for the correlation function, in terms of 
very elementary functions. This was achieved by completely diagonalizing the Fokkcr-Planck operator 
and hence writing the correlation function as a series, whose sum can be explicitly determined. 

In contrast, the motion on the three-dimensional sphere is more involved. This is largely due to the 
more complicated structure of the rotation group in three dimensions. Quite generally, the components 
of n perpendicular to a; rapidly rotate, hence decorrelate, whereas the component parallel to a; remains 
unchanged, at least when a; is constant. In the limit when the persistance angle is large, this leads to 
a two time-scale dynamics: a fast decorrelation of n is observed, corresponding to the components of 
n perpendicular to a;, followed by a much slower decorrelation of n, corresponding to the component 
parallel to ijJ. Whereas the fast decorrelation can be understood quantitatively by using elementary 
considerations, the description of the decorrelation of n at large times requires a determination of the 
largest eigenvalue of the Fokker-Planck operator. We have computed here the eigenvalues relevant to 
the long term evolution of (n(0) • n(t)). Interestingly, in the large /3 limit, the problem reduces to a 
quantum harmonic oscillator with a irrational angular momentum. 

The Ornstein-Uhlenbeck model is closely related to the equation describing the orientational degrees 
of freedom of small particles in turbulent flows , and numerical studies show that the vorticity 

of turbulent flows also has an exponential correlation [TTlfTS] . However, we have observed here that 
caution should observed in applying our model to rotation by turbulence. In two dimensions we showed 
that the correlation function of n takes a very different form when the angular velocity is generated 
by a telegraph process, despite the fact that the correlation function of a; are identical to our model. 

In summary, the notion of the 'persistence angle' introduced here appears to be the most relevant 
parameter characterising random rotation. The Ornstein-Uhlenbeck model is the simplest description 
of random motion on a sphere, and it will surely find signoficant applications, beyond the example 
considered at the end of section O However, it is a poor model for rotations of small bodies driven by 
hydrodynamic turbulence |18| . 
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7 Appendix A: Matrix representation 

7.1 Decomposition and projection. 

The aim of this appendix is to project the system of partial differential equations (|87p on the complete 
set of eigenfunctions of the spherical harmonic oscillator operator. The operators L/, introduced in 
equation ([55]) correspond to the radial part of the equation for the harmonic oscillator operators, 
equations (|681 169p associated with angular momentum quantum number Z = 0, 1 and 2 respectively. 

The eigenvalues of L/ are thus — (2n-|- Z), n being the quantum number characterizing energy, and 
the corresponding radial eigenfunction being: 



where L„ is the generalized Laguerre polynomial [7], and the normalization constant is given by 
equation (jlOTp . The following discussion uses properties of the generalised Laguerre polynomials which 
arc discussed in pS] . 

For each value of I, these eigenstates are orthogonal to each other, in the sense that: 



4^{r)=Mni exp(-rV2) r' 4'+'/'^^') 



(123) 



'nm 



(124) 
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One can therefore expand the functions m series of the 4>^: 

oo 

^(^\r) = J2^n 4\r) . (125) 

n=0 

This can then be inserted into the set of equations (|107p . Then, the equation corresponding to angular 
momentum I is projected on the set of modes (f)n\r). This leads to: 

(0(2), L^V*')) + x/V3/3(0i'),r^(^') = A(</.(2), (126) 

Clearly, 

(0W,L,V'(')) = -a(" (2n + 0. (127) 

The scalar products of {(j>^\ rip^^^ ) and {4'n \r'ip^^''), on one hand, and (0L^\r^(2)) ^^^^ {(j3'^\r^W) 
on the other hand involve the calculations of the two integrals: 

poo 

/O'l ^ / exp(-r2)r0+i x r x r24V2) (,2)^(3/2) (^2) 



and: 



= t: / cxp(-u)u5/2i(,3/2)(u)L(5/2)(y) . (129) 
2 Jo 

These two integrals can be easily computed, by using the following relation between generalized La- 
guerre polynomials: 

L(r)(x) = Lir'H^) - Lt\'\^) (130) 
and the normalization integral fj73p . One thus finds: 

I^L = 1 1 eM-uV"LT\u){L^^/^\u) - 4l/?V)) d. 
r(5/2 + m) 

= (dri.m - dn-l.m) j (lol) 

and 



1 

\ / exp(-u)u3/2i(i/2)(,)43/2)(^) (128) 



/i;^ ^ / exp(-r2)ri+2 x r x r^^ff Z^) (^2)45/2) (,2) 



71,2 

n.m 



1 Z*^ 

i exp(-.).^/24/2)(u)(4/2)(.) - 4_/f (.)) d. 

, r(7/2 + m) 
2 to! 

Using these results, the set of equations (|126p reduces to a matrix equation, with a relatively simple 
(band-) structure, as we explain below. 



23 



7.2 Matrix equations 



It is now a simple matter to rewrite the matrix equations for the quantities an\ defined in equation 



1251) . Specifically, from Eq. (|126|) . one obtains the system of equations: 

oo 

- 2na([') - 2/^3/3 A{n, m)a'^^^ = A"(°) 

OO 

2n + l)a(,i) + 2^3/3 B{n, m) a^^^ - 

(2n + 2)a(f) 



m=0 



m=0 
oo 

m=0 



(133) 



m=0 



where 



C{n, m) 



J0,1 

rt.m 



A/'„lA/'m2 



i?(ri, m) 



Z?(ri, m) 



/1,2 



(134) 



In fact, in view of the structure of the scalar products Eq. (|131ll32p . the matrix equations Eq. (|133p 
contain in fact very few terms. Explicitly, 



- 2na(J') - 2(i/Vi[(A{n,n) a^i) - A{n,n- 1) a^^\) 
-(2n+l)ali) +/3/V3 [2(5(71, n) a'^^^ ~ B{n,n + I) a|°|i 
-{C{n,n) a(f)-C(n,n-l)a|f2i) 
^(2n + 2)a,[f) + ^273/3 [(^'(n, n) a^^) - D{n, n + 1) a^^ji) 



= Aai^) 



(135) 



The structure of the matrix defined by equations (jl35p . is easy to program in a routine that can 
diagonalize a real matrix, with the infinite-dimensional matrix truncated to a finite size by including 
only coefficients with n < nmax- Although the structure of the matrix is very sparse - the matrix 
a simple band structure - a general routine from NAG has been used to determine the spectrum. We 
checked that the eigenvalues converge as the number of coefficients rimax increases. We were able to 
obtain the eigenvalues with the largest real parts, i.e., the one that correspond to the slowest decay of 
the correlation function. 



8 Appendix B: Telegraph-noise model 

When solving physical problems it is often tacitly assumed that the correlation function of a stochastic 
signal is sufficient to characterise its properties. For example, if we were to use a different stochastic 
process to generate the angluar velocity, we might expect that the correlation function C{t) = (n(i) • 
n(0)) would be little changed if the correlation function of the angular velocity reamins the same, 
namely {ijj(t)u}{<d)) = (D/7) exp(— 7|t|). This assumption may not be valid in general. To illustrate 
this point, here we determine the correlation function ^ for random motion on a circle in the case 
where the angular velocity is generated by a telegraph noise process, which has precisely the same 
exponential correlation function as the Ornstein-Uhlenbeck process. The correlation function of the 
direction vector, defined by ([H]), is found to be very different from ([55]). 

In the telegraph noise model the angular velocity Lo{t) takes just two discrete values, ibwo (where 
Wo is a constant). The angular momentum makes random transitions between these values, with a rate 
constant R (so that the probability of transition in a short time interval of length 5t is 5P = RSt). 
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The rotation angle 9 satisfies 9 = Lo{t), as before. Define P±{9,t) to be the probability density to be 
located at 9 at time t, with uj{t) = iwo- The probability densities satisfy 



dP± 



± Wo 



dP± 



RP^ - RP± 



(136) 



Formally, this equation can be written dt\P) = -^l^*): where \P) is a function vector representing 
{P+{9,t), P-{9,t)). This equation can be solved by seeking eigenfunctions of in the form = 
exp(in6')(c+, c_). The vector c = (c+,c_) is an eigenvector of the 2x2 matrix 



-inujQ — R R 
R incjQ ^ R 



This matrix has eigenvalues 



All 



\n± — — _i_ y j^L — (I, ujQ 

The general solution of p36p can be expressed as a linear combination of eigenfunctions: 



(137) 
(138) 



P+{9,t) 
P-{0,t) 



n— — oo 



i? 



+ a„_ exp(A„_t + i7i6') 



inwo + ^ R? — n^cijQ 
i? 



In order to evaluate the correlation function (|6]) we must compute C(t) = (cos( 
condition P±{9) = S{9)/2. Thus 



(139) 

with the initial 



C{t) 



d9 cos9[P+{9,t) + P^{9,t)] 



(140) 



= 27rRe 



ai+ [ R + iujQ + J R^ - w§ exp(Ai+i) + ai_ i? + iwo - J R'^ - cuq exp(Ai_t) 



The coefficients a„± are easily determined from the initial distribution: 



an± 



v/i?2 _ „2^2 ± _ i^^g) 



SttR^R^ - 



,2,, ,2 



and hence 



Cit) = -[cxp(A+t)+cxp(A_f)] + 



i? 



:[exp(A+t) - exp(A_t)] 



(141) 



(142) 



where X± = Xi± = — i?± R"^ ~ ^o- similar and somewhat simpler calculation gives the correlation 
function of Lj{t): 

(w(t)a;(0)) = w^exp(-2i?|i|) . (143) 

Because this correlation function has the same structure as that of the Ornstein-Uhlenbeck process, 
we can define the parameters luq, R of the telegraph noise model in terms of the parameters 7, D of 
the Ornstein-Uhlenbeck process: by comparison of ()143p with Q we have 



7 - 2i? , 2Rwl , /? = ^ . 



(144) 



Expressed in terms of the same variables as the Ornstein-Uhlenbeck process, the correlation function 
of the telegraph noise model is 



C(t) = - — , exp 
2^1-4/32 



1- yr^^ 



It 



yi^4^ 



1 



2v/r^4^ 



■ exp 



1 + v/1 - 4/32 



It 



(145) 



This correlation function is significantly different from ([55]): for example (|145p is oscillatory when 
(3 > 1/2. 
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